Filogenetska stabla

Filogenija proučava evolutivni razvoj organizama. Za reprezentaciju evolutivnog razvitka vrsta i njihovih udaljenosti od zajedničkih predaka koriste se filogenetska stabla. Za konstrukciju takvog stabla, potrebno je definisati meru koja će dobro oslikavati udaljenosti između vrsta. Jedna takva mera može biti broj različitih nukleotida na istim pozicijama homologih fragmenata DNK/RNK sekvenci ili aminokiselina u proteinskim sekvencama. Konstrukcija stabala kao preduslov ima upotrebu aditivne mere za meru udaljenosti između elemenata. Aditivna mera udaljenosti je ona mera koja obezbeđuje da zbir svih grana na putu između dva čvora stabla odgovara udaljenosti polaznog i krajnjeg čvora prema zadatoj meri. Iako evolucija nalaže da je stablo filogenije korensko stablo (postoji jedan zajednički predak (npr. vrsta) da sve žive organizme na planeti), u zavisnosti od analize, filogenetsko stablo ne mora biti korensko.

Čvorovi stabla mogu biti unutrašnji i spoljašnji. Spoljašnji čvorovi, listovi, odgovaraju uzorcima između kojih se analizira filogenija (vrstama, ćelijama, proteinima, ...). Unutrašnji čvorovi predstavljaju hipotetičke zajedničke pretke uzoraka čiji je roditelj u stablu unutrašnji čvor.

Dužine grana stabla predstavljaju udaljenost vrsta koja je direktno proporcionalna vremenskom periodu koji je protekao od razdvajanja uzoraka od zajedničkog pretka. Grana koja spaja list i unutrašnji čvor zvaće se limb. Kako udaljenost unutrašnjeg čvora od ostalih čvorova nije poznata, dužina ove grane može se izračunati korišćenjem jednog čvora čiji je navedeni unutrašnji čvor takođe sused i jednog čvora čiji navedeni unutrašnji čvor nije sused:

Čvorovi $j$ i $i$ imaju zajedničkog suseta $m$, dok čvor $k$ nije sused unutrašnjeg čvora $m$. Udaljenost čvora $j$ od svog unutrašnjeg suseda $m$ jednaka je $D_{j,m} = (D_{j,k} * D_{i,k} - D_{i,j}) / 2$ (slika iznad).

Unutrašnji čvor nije poznat a potrebno je na neki način utvrditi koji je to čvor $i$ koji ima zajedničkog unutrašnjeg suseda sa čvorom $j$. Traženi čvor $i$ je onaj čvor koji daje minimum izraza $D_{j,m} = (D_{j,k} * D_{i,k} - D_{i,j}) / 2$, jer se time garantuje da se $i$ nalazi na najmanjoj udaljenosti od čvora $j$ u odnosu na unutrašnji čvor $m$ i udaljeni čvor $k$. Iz navedenog izraza, lako se mogu izračunati i ostale udaljenosti između navedenih čvorova.

Neka je aditivna $D$ matrica udaljenosti između listova:

In [1]:
D = [[ 0, 13, 21, 22],
     [13,  0, 12, 13],
     [21, 12,  0, 13],
     [22, 13, 13,  0]]

stabla:

Dužina limb grane čvora $j$ računa se pomoću funkcije (argument $n$ predstavlja ukupan broj listova):

In [2]:
def limb(D, j, n):
    
    min_i = None
    min_k = None
    
    min_length = float('inf')
    
    for i in range(n):
        for k in range(n):
            if i != j and k != j:
                length = (D[j][k] + D[i][j] - D[i][k]) / 2
                if length < min_length:
                    min_length = length
                    min_i = i
                    min_k = k
                    
    return (min_i, min_k, min_length)
In [3]:
i = 0
j = 1
n = len(D)
Dim = limb(D, i, n)
Djm = limb(D, j, n)
print(f'i = {i}, j = {j}')
print(f'D(i, m) = {Dim[2]}, preko čvorova j = {Dim[0]} i k = {Dim[1]}')
print(f'D(j, m) = {Djm[2]}, preko čvorova j = {Djm[0]} i k = {Djm[1]}')
i = 0, j = 1
D(i, m) = 11.0, preko čvorova j = 1 i k = 2
D(j, m) = 2.0, preko čvorova j = 0 i k = 2

Filogenetsko stablo će biti predstavljeno kao usmereni graf:

In [4]:
import copy

class Graph:
    # Inicijalizacija listom susedstva
    def __init__(self, adjacency_list):
        self.G = adjacency_list

    # Funkcija vraća listu susedstva grafa 
    def adjacency_list(self):
        return copy.deepcopy(self.G)
        
    # Dodavanje novog cvora
    def add_vertex(self, v):
        if v not in self.G:
            self.G[v] = []
        
    # Lista suseda čvora v
    def get_neighbors(self, v):
        if v in self.G:
            return self.G[v]
        else:
            []
            
    # Dodavanje novog suseda zadatog čvora na zadatoj udaljenosti
    def add_neighbor(self, vertex, neighbor, distance):
        self.G[vertex].append((neighbor, distance))

        if neighbor not in self.G:
            self.G[neighbor] = []

        self.G[neighbor].append((vertex, distance))
        
    # Uklanjanje susednog čvora i grane do istog
    def remove_neighbor(self, vertex, neighbor, distance):
        self.G[vertex].remove((neighbor, distance))
        self.G[neighbor].remove((vertex, distance))
            
    # Pronalaženje puta između dva čvora (iterativna DFS)
    def find_path(self, source, destination):
        stack = [source]
        visited = set([source])

        while len(stack) > 0:
            v = stack[-1]

            if v == destination:
                return stack

            has_neighbors = False

            for (w, weight) in self.get_neighbors(v):
                if w not in visited:
                    stack.append(w)
                    visited.add(w)
                    has_neighbors = True
                    break

            if not has_neighbors:
                stack.pop()

        print('Path not found')
        return []
    
    # Udaljenost između susednih čvorova
    def distance_between_neighbors(self, vertex_i, vertex_j):
        neighbors = self.get_neighbors(vertex_i)

        for (w, weight) in neighbors:
            if w == vertex_j:
                return weight

        return None
    
    # Dodavanje čvora na putu između dva čvora, na zadatoj udaljenosti od polaznog čvora
    def add_vertex_on_path(self, source, destination, distance):
        path = self.find_path(source, destination)

        i = 0
        j = 1

        vertex_i = path[i]
        vertex_j = path[j]

        current_distance = self.distance_between_neighbors(vertex_i, vertex_j)

        while current_distance < distance:
            i += 1
            j += 1

            vertex_i = path[i]
            vertex_j = path[j]

            current_distance += self.distance_between_neighbors(vertex_i, vertex_j)

        if current_distance == distance:
            return vertex_j

        else:
            distance_between = current_distance - distance
            return self.new_vertex_between_vertices(vertex_j, vertex_i, distance_between)
        
        
    # Dodavanje čvora između dva susedna čvora na zadatoj udaljenosti od polaznog čvora
    def new_vertex_between_vertices(self, vertex_i, vertex_j, distance_between):
        new_vertex = 'X{}{}'.format(vertex_i, vertex_j)

        distance_ij = self.distance_between_neighbors(vertex_i, vertex_j)
        distance_diff = distance_ij - distance_between

        self.remove_neighbor(vertex_i, vertex_j, distance_ij)

        self.add_neighbor(vertex_i, new_vertex, distance_between)
        self.add_neighbor(new_vertex, vertex_j, distance_diff)

        return new_vertex
In [5]:
# PRIMER GRAFA
adjacency_list = {
    'a': [('b', 3)],
    'b': [('a', 3), ('c', 2), ('d', 1)],
    'c': [('b', 2)],
    'd': [('b', 1)]
}

G = Graph(adjacency_list)

print("get_neighbors('a'): ", G.get_neighbors('a'))
print("find_path('a', 'c'): ", G.find_path('a', 'c'))
print("distance_between_neighbors('a', 'b'): ",G.distance_between_neighbors('a', 'b'))

new_node = G.new_vertex_between_vertices('b', 'c', 1.5)
print("Novi čvor: ",new_node)
print("adjacency_list(): ",G.adjacency_list())
get_neighbors('a'):  [('b', 3)]
find_path('a', 'c'):  ['a', 'b', 'c']
distance_between_neighbors('a', 'b'):  3
Novi čvor:  Xbc
adjacency_list():  {'a': [('b', 3)], 'b': [('a', 3), ('d', 1), ('Xbc', 1.5)], 'c': [('Xbc', 0.5)], 'd': [('b', 1)], 'Xbc': [('b', 1.5), ('c', 0.5)]}

Konstrukcija stabla aditivne filogenije pomoću strukture grafa:

In [6]:
def additive_phylogeny(D, n):
    # Ako stablo ima dva čvora, oni čine granu stabla
    # dok je njihova međusobna udaljenost navedena u matrici udaljenosti
    if n == 2:
        return Graph({
            0: [(1, D[0][1])],
            1: [(0, D[0][1])]
        })
    
    # U suprotnom, odabira se proizvoljni čvor j (u ovom slučaju n)
    # i određuje se njegova udaljenost do svog unutrašnjeg suseda (limb_length)
    (i, k, limb_length) = limb(D, n - 1, n)
    
    # Zatim se preračunava matrica udaljenosti tako da oslikava udaljenosti u stablu
    # bez lista n (od svih udaljenosti između čvora n i ostalih čvorova 
    # oduzima se udaljenost čvora n od svog unutrašnjeg čvora)
    #
    # Neformalno: Čvor n se pomera tako da se preklapa sa svojim unutrašnjim susedom
    # što daje stablo u kome je čvor n uklonjen
    for j in range(n - 1):
        D[j][n - 1] -= limb_length
        D[n - 1][j] = D[j][n - 1]
        
    # Din predstavlja udaljenost čvora i od sada zajedničkog unutrašnjeg čvora
    Din = D[i][n - 1]
    
    # Rekurzivno se konstruiše stablo bez čvora n...
    T = additive_phylogeny(D, n - 1)
    
    # ... a zatim se u stablo dodaju unutrašnji čvor m, na udaljenosti Din,
    # na putu između čvorova i, k...
    m = T.add_vertex_on_path(i, k, Din)
    
    # i konačno, prethodno izvačeni čvor n kao sused unutrašnjeg čvora m na udaljenosti limb_length
    T.add_neighbor(m, n-1, limb_length)
    
    return T
In [7]:
D = [[ 0, 13, 21, 22],
     [13,  0, 12, 13],
     [21, 12,  0, 13],
     [22, 13, 13,  0]]

# Broj čvorova
n = 4

T = additive_phylogeny(D, n)
print(T.adjacency_list())
{0: [('X10', 11.0)], 1: [('X10', 2.0)], 'X10': [(1, 2.0), (0, 11.0), ('X2X10', 4.0)], 2: [('X2X10', 6.0)], 'X2X10': [(2, 6.0), ('X10', 4.0), (3, 7.0)], 3: [('X2X10', 7.0)]}
In [8]:
# OPCIONO: Vizualizacija korišćenjem networkx biblioteke (pip3 install networkx)
import networkx as nx
import matplotlib.pyplot as plt

def graph_to_nx(G):
    Gx = nx.Graph()
    Gx.add_nodes_from(list(G.keys()))

    for v in G:
        for w in G[v]:
            Gx.add_edge(v, w[0], weight= w[1])

    nx.draw(Gx, with_labels=True, font_weight='bold')

    plt.show()
    
    return Gx

Tx = graph_to_nx(T.adjacency_list())

    
<Figure size 640x480 with 1 Axes>

Drugi pristup konstrukciji filogenetskog stabla je UPGMA algoritam. Inicijalno, svi elementi (listovi stabla) predstavljaju klastere sa po jednim elementom. Iterativno se klasteri na najmanjim međusobnim udaljenostima spajaju u nove klastere, sve dok se ne dođe do jednog klastera koji sadrži sve elemente. Svako spajanje klastera predstavlja unutrašnji čvor filogenetskog stabla. Udaljenosti između klastera računaju se kao srednja vrednost međusobnih rastojanja svih elemenata iz jednog klastera od elemenata drugog klastera. Za razliku od aditivne filogenije, UPGMA algoritam konstruiše korensko stablo.

Struktura klastera:

In [9]:
class Cluster:
    def __init__(self, elements):
        self.elements = elements
        self.age = 0
        self.left = None
        self.right = None
        self.label = str(elements)
        
    def add_children(self, left, right):
        self.left = left
        self.right = right
        
    # Nivo klastera u hijerarhiji
    def update_age(self, age):
        self.age = age
        
    # Udaljenost klastera od drugog klastera
    def distance(self, D, cluster_2):
            d = 0
    
            n = len(self.elements)
            m = len(cluster_2.elements)

            for i in self.elements:
                for j in cluster_2.elements:
                    d += D[i][j]

            d = d / (n * m)
            return d
        
    # Formatiranje klastera za ispis
    def __str__(self):
        return 'Label: {}, Age: {}'.format(self.label, self.age)

Nalaženje dva klastera, u nizu klastera, na najmanjoj udaljenosti:

In [10]:
def two_closest(clusters, D):
    min_i = None
    min_j = None
    
    min_distance = float('inf')
    
    for cluster_i in clusters:
        for cluster_j in clusters:
            if cluster_i != cluster_j:
                current_distance = cluster_i.distance(D, cluster_j)
                if current_distance < min_distance:
                    min_distance = current_distance
                    min_i = cluster_i
                    min_j = cluster_j
                    
    return (min_i, min_j, min_distance)

Transformacija listova stabla u klastere

In [11]:
def leaves_to_clusters(D):
    n = len(D)
    clusters = []
    
    for i in range(n):
        clusters.append(Cluster([i]))
        
    return clusters
In [12]:
# PRIMER
clusters = leaves_to_clusters(D)
for c in clusters:
    print(c)

print()
print(D)
print()
(c1, c2, dist) = two_closest(clusters, D)
print('C1: ')
print(c1)
print('C2: ')
print(c2)
print('Distance: {}'.format(dist))
Label: [0], Age: 0
Label: [1], Age: 0
Label: [2], Age: 0
Label: [3], Age: 0

[[0, 13, 11.0, 15.0], [13, 0, 2.0, 6.0], [11.0, 2.0, 0, 6.0], [15.0, 6.0, 6.0, 0]]

C1: 
Label: [1], Age: 0
C2: 
Label: [2], Age: 0
Distance: 2.0
In [13]:
def UPGMA(D, n):
    
    # Inicijalizacija jednočlanih klastera, listova stabla
    clusters = leaves_to_clusters(D)
    T = Graph({})
    for c in clusters:
        T.add_vertex(c.label)
        
    while len(clusters) > 1:
        # Pronalaženje najbližih klastera
        (c_i, c_j, dist) = two_closest(clusters, D)
        
        # Spajanje najbližijh klastera
        c_new = Cluster(c_i.elements + c_j.elements)
        c_new.left = c_i
        c_new.right = c_j
        c_new.update_age(dist / 2)
        
        # Dodavanje novog čvora u stablo
        T.add_vertex(c_new.label)
        T.add_neighbor(c_new.label, c_i.label, dist/2)
        T.add_neighbor(c_new.label, c_j.label, dist/2)
        
        new_clusters = [c_new]
        
        # Ažuriranje skupa klastera
        for c in clusters:
            if c != c_i and c != c_j:
                new_clusters.append(c)
                
        clusters = new_clusters
        
    root = clusters[0]

    return (T, root)
In [14]:
# Primer

D = [[ 0, 13, 21, 22],
     [13,  0, 12, 13],
     [21, 12,  0, 13],
     [22, 13, 13,  0]]

n = 4

graph, root = UPGMA(D, n)

# print(graph.adjacency_list())
_ = graph_to_nx(graph.adjacency_list())

Koren dobijenog stabla je čvor sa oznakom svih elemenata, u primeru to je "[1,2,3,0]"

In [15]:
print(root.age)
9.333333333333334

Još jedan aloritam za konstrukciju Neighbor Joining koji spaja elemente na najmanjoj udaljenosti na osnovu $D^{*}$ matrice, čija konstrukcija se za potrebe ovog kursa može uzeti zdravo za gotovo.

Matrica D će, zbog pogodnijeg oblika, biti transformisana u mapu:

In [16]:
import copy

class D_Map:
    # Inicijalizacija nove mape
    def __init__(self, D):
        self.D = D
        n = len(D)
        self.n = n

        D_map = {}
        for i in range(n):
            D_map[i] = {}

            for j in range(n):
                D_map[i][j] = D[i][j]
                
        self.D_map = D_map

    # Formatiranje objekta za ispis
    def __str__(self):
        return f'{self.D_map}'
    
    # Lista (kljuc,vrednost) parova elemenata
    def items(self):
        return self.D_map.items()
    
    # Funkcija vraća udaljenost između čvorova (i, j)
    def dist(self, i, j):
        return self.D_map[i][j]
    
    # Dodavanje novog cvora m u odnosu na najbliži par suseda 
    # (funkcija za nalazenje najbližih suseda je u nastavku klase)
    def add_node(self, m, i, j):
        nodes = [x[0] for x in list(self.D_map.items())]
        self.D_map[m] = {}

        for k in nodes:
            if k != i and k != j:
                self.D_map[k][m] = 1/2 * (self.D_map[k][i] + self.D_map[k][j] - self.D_map[i][j])
                self.D_map[m][k] = self.D_map[k][m]

        self.D_map[m][m] = 0
        self.n += 1
        
    # Uklanjanje cvora iz mape
    def remove_node(self, removed_node):
        D_map = self.D_map
        del D_map[removed_node]

        for node in D_map:
            if removed_node in D_map[node]:
                del D_map[node][removed_node]
                
        self.n -= 1
                
    # Zbir svih udaljenosti čvora do ostalih čvorova
    def total_distance(self, node_i):
        total = 0

        for node_j in self.D_map[node_i]:
            total += self.D_map[node_i][node_j]

        return total
    
    # Konstrukcija D* mape
    def map_star(self):
        D_map = self.D_map
        D_map_star = copy.deepcopy(D_map)
        n = self.n

        for i in D_map:
            for j in D_map[i]:
                if i != j:
                    D_map_star[i][j] = (n - 2) * D_map[i][j] - self.total_distance(i) - self.total_distance(j)
                else:
                    D_map_star[i][j] = 0

        return D_map_star
    
    # Pronalaženje najbliža dva elementa u odnosu na D*
    def minimal_pair(self):
        min_distance = float('inf')
        min_i = None
        min_j = None
        
        D_star = self.map_star()

        for i in D_star:
            for j in D_star[i]:
                if i != j:
                    current_distance = D_star[i][j]
                    if current_distance < min_distance:
                        min_i = i
                        min_j = j
                        min_distance = current_distance

        return (min_i, min_j)
In [17]:
# Primer:
D = [[ 0, 13, 21, 22],
     [13,  0, 12, 13],
     [21, 12,  0, 13],
     [22, 13, 13,  0]]

dmap = D_Map(D)
print(dmap)
{0: {0: 0, 1: 13, 2: 21, 3: 22}, 1: {0: 13, 1: 0, 2: 12, 3: 13}, 2: {0: 21, 1: 12, 2: 0, 3: 13}, 3: {0: 22, 1: 13, 2: 13, 3: 0}}
In [18]:
def neighbor_joining(D):
    if D.n == 2:
        (i, j) = [x[0] for x in list(D.items())]
        return Graph({
            i: [(j, D.dist(i,j))],
            j: [(i, D.dist(j,i))]
        })

    (i, j) = D.minimal_pair()

    delta = (D.total_distance(i) - D.total_distance(j)) / (n - 2)
    limb_length_i = 1/2 * (D.dist(i,j) + delta)
    limb_length_j = 1/2 * (D.dist(i,j) - delta)

    m = '{}+{}'.format(i, j)

    D.add_node(m, i, j)
    D.remove_node(i)
    D.remove_node(j)

    T = neighbor_joining(D)

    T.add_vertex(m)
    T.add_neighbor(m, i, limb_length_i)
    T.add_neighbor(m, j, limb_length_j)

    return T
In [19]:
# Primer:
D = [[ 0, 13, 21, 22],
     [13,  0, 12, 13],
     [21, 12,  0, 13],
     [22, 13, 13,  0]]

D_map = D_Map(D)
T = neighbor_joining(D_map)
print(T.adjacency_list())
_ = graph_to_nx(T.adjacency_list())
{'0+1': [('2+3', 4.0), (0, 11.0), (1, 2.0)], '2+3': [('0+1', 4.0), (2, 6.25), (3, 6.75)], 2: [('2+3', 6.25)], 3: [('2+3', 6.75)], 0: [('0+1', 11.0)], 1: [('0+1', 2.0)]}